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Using ah initio calculations, we study the electronic and structural properties of vacancies and 
hydrogen adsorbates on trilayer graphene. Those defects are found to share similar low-energy 
electronic features, since they both remove a pz electron from the honeycomb lattice and induce a 
defect level near the Fermi energy. However, a vacancy also leaves unpaired a electrons on the lattice, 
which lead to important structural differences and also contribute to magnetism. We explore both 
ABA and ABC stackings and compare properties such as formation energies, magnetic moments, 
spin density and the local density of states (LDOS) of the defect levels. These properties show a 
strong sensitivity to the layer in which the defect is placed and smaller sensitivities to sublattice 
placing and stacking type. Finally, for the ABC trilayer, we also study how these states behave in the 
presence of an external field, which opens a tunable gap in the band structure of the non-defective 
system. The pz defect states show a strong hybridization with band states as the field increases, 
with reduction and eventually loss of magnetization, and a non-magnetic, midgap-like state is found 
when the defect is at the middle layer. 


I. INTRODUCTION 

Graphene, a 2D material composed of carbon atoms 
arranged in a honeycomb lattice, has attracted much at¬ 
tention from the scientific community in the last decade 
due to its unique electronic and structural properties and 
its potential applications in a wide variety of fields mm- 
Among many topics studied in graphene research, one 
of particular interest is the study of defects, and how 
they affect the electronic and structural properties of the 
sheet. In particular, two well-studied types of defects in 
graphene are vacancies and hydrogen adsorbates on top 
of the sheet, which can be produced experimentally by 
bombarding the sheet with high-energy particles such as 
protons mm- Such defects generally induce a local mag¬ 
netism in the graphene sheet, a highly desirable prop¬ 
erty for spintronics applications mm- Magnetic materi¬ 
als commonly used in current technological applications 
are based on heavy elements, with d and / orbitals, so 
the realization of s and p electron magnetism in graphene 
is very attractive. 

In single layer graphene, isolated vacancies and hydro¬ 
gen adsorption are found to have similar low-energy prop¬ 
erties. Both defects remove a Pz electron from the hon¬ 
eycomb lattice, inducing a ” quasi-localized” defect level 
near the Fermi energy. The associated wavefunction is 
found to be long-ranged, decaying as 1/r, where r is 
the distance to the defect center, as predicted by the¬ 
ory [HI [9] and revealed by STM experiments on vacancies 
in graphite uni. 

Pure trilayer graphene also attracts increasing atten¬ 
tion due to its stacking-dependent electronic proper¬ 
ties. The two main types of stacking, ABA and ABG, 
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are shown in Fig. In ABA (or Bernal) stacking, 
the top and bottom layers have the same orientation, 
while the middle layer is rotated by 60° with respect to 
them. The associated low-energy band structure shows 
a pair of single-layer and bilayer-like bands, with linear 
and quadratic dispersions, respectively mHH]. On the 
other hand, in ABG stacking, all the layers are rotated 
by the same angle with respect to each other, so they 
all have different orientations. The corresponding low- 
energy band structure in this case shows a pair of bands 
with cubic dispersion. In addition, the application of 
an external electrical field perpendicular to the sheets 
in ABG trilayer breaks its inversion symmetry and in¬ 
duces a tunable band gap, a very desirable property for 
electronic devices [m na HE]. Point defects in trilayer 
graphene are expected to give rise to a wider variety of 
electronic and structural features, since they might de¬ 
pend on the stacking order and the site where the defect 
is. 

In this work, we explore these different possibilities, by 
employing ab-initio calculations based on Density Func¬ 
tional Theory (DFT). We compare the features of the 
low-energy defect level induced by an isolated vacancy 
and by adsorption of a single hydrogen atom on top of 
ABA and ABG trilayer graphene. Specifically, we place 
the defect at different non-equivalent sites and study how 
properties such as formation energies, magnetizations, 
band structure, spin densities and local density of states 
(LDOS) behave in each case. For the ABG trilayer, we 
also study the effect of an external electrical field, where 
we look in particular for the possibility of a mid-gap de¬ 
fect level, as already pointed out by a previous tight- 
binding calculation m- Our paper is organized as fol¬ 
lows: in the next Section, we present the methodology 
used in this work. In Sections [ni| and [IV| we discuss our 
results for hydrogen adsorption and isolated vacancies, 
respectively, where the properties mentioned above are 



2 


ABA 




FIG. 1: Top view of ABA (top) and ABC (bottom) trilayer 
graphene. Yellow and dark yellow spheres represent carbon 
atoms on top and middle layers, respectively. The basis atoms 
are labeled by their sublattice (A or B) and layer indexes (1, 
2 or 3). Sites with more than one index indicate different 
carbon atoms that are on top of each other. The bonds in the 
top, middle and bottom layers are represented by the black, 
red and blue lines, respectively (blue is only visible in the 
ABC trilayer). 

discussed in different subsections. Finally, in Section [V| 
we summarize our results and present our conclusions. 


II. METHODOLOGY 

We use first-principles calculations based on the Den¬ 
sity Functional Theory (DFT) (TSUTOj, as implemented in 
the SIESTA code [20] . We use norm-conserving Troullier- 
Martins pseudopotentials [21] for the ion-electron inter¬ 
action, a PBE-GGA exchange-correlation functional [22] 
and a double-C-polarized (DZP) pseudoatomic basis set 
for the expansion of the electronic wavefunctions. The 
real space grid energy cutoff is set to 150 Ry and we use 
a 6 X 6 X 1 Monkhorst-Pack k-point grid [23] and an elec¬ 
tronic smearing (Eermi-Dirac-like function) of 100 K in 
all calculations. 

Eor all configurations considered, we let the atomic co¬ 
ordinates relax in all directions until all the forces are 
smaller than 0.04 eV/A. However, for the cases where an 


electrical field is applied (ABC trilayer only), we use the 
optimized coordinates from the corresponding configura¬ 
tion in zero field. In order to simulate the vacancy and 
the hydrogen adsorption defects, we use supercells of di¬ 
mension 6a X 6a, where a = 2.47 A is the lattice constant 
of pure graphene. This corresponds to a defect density 
of ^ 5 X 10^^ cm“^, which is a high density compared 
to current experimental values m- The optimized value 
for the interlayer distance is 3.27 A in all cases and we 
use a vacuum of 10 A in the direction perpendicular to 
the layers in order to avoid interactions between periodic 
images [35] . We have tested supercells of larger sizes for a 
few configurations and we didn’t see qualitative changes 
of our results. The different configurations studied in this 
work are labeled by the stacking type (ABA or ABC), the 
type of defect and the site where the defect is, following 
the notation of Eig. Eor hydrogen adsorption, the H 
atom always lies on top of the nearest carbon atom and 
the bond length is relaxed together with the the struc¬ 
ture. Eor A1 and B1 sites, we always choose the hydrogen 
atom to be outside, instead of between the layers as in 
the A2 and B2 sites. We performed both spin polarized 
and unpolarized calculations, but we report mainly our 
results for the spin polarized case, since they are of more 
interest, unless mentioned otherwise. 


III. HYDROGEN ADSORPTION 

A. Structural properties, formation energies and 
magnetic moments 

Our calculations show that the structural changes due 
to hydrogen adsorption are small and nearly stacking- 
independent. The hydrogen atom ’’pulls” the carbon 
atom it binds to from the layer, along with some of its 
neighbors. This effect is due to a local sp^ — sp^ rehy¬ 
bridization in order to form the C-H bond. As shown 
in Table |T] the relaxed length of the C-H bond lies be¬ 
tween 1.13 — 1.15 A in all cases. However, the magnitude 
of the displacement of the carbon atoms depends on the 
defect site. When the H atom is on top of A1 or B1 
sites (outside the trilayer), the out-of-plane displacement 
of the C atom in the C-H bond ranges from 0.41 to 0.49 
A, while its nearest neighbors move between 0.07 — 0.14 
A. On the other hand, when the H atom is on top of 
A2 and B2 sites, between the top and middle layers, the 
displacements are smaller: between 0.20 — 0.28 A and 
—0.02 — 0.00 A, respectively. We believe this effect is re¬ 
lated to the interlayer interaction: in the first case, the 
carbon atoms have more freedom to relax, while in the 
second case they are tightly squeezed between the two 
adjacent layers, leading to a smaller displacement. We 
also do not see any changes in the interlayer distances 
due to the presence of the hydrogen atom. Einally, we do 
not see any noticeable structural changes in the adjacent 
pristine layers. 

In order to further study the sensitivity to the defect 
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TABLE I: Results for spin-polarized calculations of hydrogen adsorption on trilayer graphene: C-H bond length (dc-rr), 
out-of-plane displacements of the carbon atom in the C-H bond (dcn) its first neighbors (dci), formation energies (Ef) 
and magnetic moments per unit supercell (M) for each configuration studied. For dci , negative values represent displacements 
in a direction opposite to the H atom. 


Defect Site 

dc-H (A) 

dcH (A) 

dci (A) 

Ef (eV) 

M {hb) 

ABA trilayer 

Al 

1.13 

0.41 

0.07 

1.14 

1.00 

A2 

1.15 

0.20 

-0.01 

1.37 

0.00 

B1 

1.13 

0.45 

0.11 

1.16 

1.00 

B2 

1.13 

0.25 

-0.01 

1.23 

0.01 

ABC trilayer 

Al 

1.13 

0.48 

0.13 

1.11 

0.91 

A2 

1.14 

0.24 

-0.02 

1.31 

0.51 

B1 

1.14 

0.49 

0.14 

1.17 

1.04 

B2 

1.13 

0.28 

0.00 

1.29 

0.87 


sites, we calculate the formation energies and magnetic 
moments for each configuration. For hydrogen adsorp¬ 
tion, the formation energy is calculated as 

Eh2 

^trilayer-\-H ^trilayer ^ 5 ( 1 ) 

where Etriiayer+H IS the total energy of a given configu¬ 
ration, Etriiayer IS the total energy of the corresponding 
pristine trilayer (with the same supercell size) and Eh 2 
is the total energy of an isolated hydrogen molecule, cal¬ 
culated using a large cubic supercell. The results are 
also reported in Table [Ij As expected, the formation en¬ 
ergies for H adsorption in middle layer sites are larger 
than those for the external layers. The highest value is 
found for the A2 site in the ABA trilayer, since this site 
has two carbon atoms on top of it, one in each external 
layer. In the ABC trilayer, the A2 and B2 sites have 
similar environments, and the energies differ slightly by 
the placement of the H atom: between two C atoms in 
the former and between a C atom and a hexagon center 
in the latter. Finally, for both trilayers, the formation 
energies for the B1 site are larger than those of the Al 
site, again due to the presence of a C atom on the top 
of the B1 site. For calculations without spin polariza¬ 
tion, the formation energies are systematically larger by 
a few tenths of meV, and follow the same trends of Ta¬ 
ble [Tj The structural parameters, formation energies and 
magnetizations reported in Table |l] are qualitatively con¬ 
sistent with previous calculations on single and bilayer 
graphene [5l[25H27]. On the other hand, a sensitivity of 
these properties to the defect placement is not seen in 
single layer graphene (due to the equivalence of the A 
and B sublattices) and there is no layer sensitivity in the 
bilayer (due to its inversion symmetry). 

The magnetization shows the most striking differences 
between the trilayers. For the ABA trilayer, when the 
defect is on one of the middle layer sites, there is no net 
magnetization in the supercell after the structural relax¬ 
ation. On the other hand, for the external layer sites, a 
’’full” magnetization is observed, coming from a full occu¬ 
pation of one of the spin-split defect levels. For the ABC 


trilayers, the Al, A2 and B2 sites show partial magneti¬ 
zation, while the B2 site retains full magnetization. We 
now discuss these results in terms of the band structure. 


B. Band structure 

The band structures for each defect configuration are 
shown in Fig. We also included the band structure 
for the pristine ABA and ABC trilayers, labeled ’’pure” 
in the corresponding frames. In a 6 x 6 supercell, the K 
and K' points from the primitive Brillouin Zone (BZ) of 
pure trilayer graphene are both folded into the F point 
of the super cell BZ, hence many double degeneracies are 
observed in the bandstructure in these cases. The intro¬ 
duction of a defect in the supercell lifts most of these 
degeneracies and also introduces a new defect level near 
the Fermi energy, which has a small bandwidth and is 
mostly flat throughout the BZ. The high density of states 
associated with such a level leads to a strong exchange 
interaction and a spin-split occurs, the so called Stoner 
effect, as is observed in most cases. We now discuss each 
case in more detail. 

In the ABA trilayer, the band structure near the Fermi 
energy is composed of a pair of bands with linear dis¬ 
persion and two pairs of bands with quadratic disper¬ 
sion, similar to those observed in single layer and bilayer 
graphene, respectively. Due to the lack of inversion sym¬ 
metry in this trilayer (it only has mirror symmetry), both 
type of bands have small gaps of a few tenths of meV. 
With the adsorption, we observe that the gap between 
the linear bands is increased to about 0.1 eV in all cases 
and the linear behavior is lost in a greater vicinity of the 
r point. On the other hand, the gap between the closest 
quadratic bands is preserved. 

In the A2 and B2 cases, where the defect is bound to 
the middle sheet, the band structure does not show spin 
polarization, while in the Al and B1 cases a spin-split 
occurs, which is specially large in the B1 case. This be¬ 
havior can be understood by looking at the bandwidth 
of the defect level in a corresponding non-polarized cal- 
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FIG. 2: Band structure results for hydrogen adsorption on ABA (top) and ABC (bottom) trilayer graphene. The H atom is 
on top of the site indicated in each frame. Spin up (down) levels are represented by black (red dashed) lines and the Fermi 
level is set to zero in all cases. Energies are in eV. For the pure trilayers (hrst column), the calculation is performed without 
spin polarization, so only black lines are shown. 


culation. For the Al and B1 sites (Fig. [^, our cal¬ 
culations show that the induced defect bands are much 
flatter than those induced in the A2 and B2 cases shown 
in Fig. Therefore, the density of states at the Fermi 
energy in these cases is larger and the Stoner criterion 
for exchange splitting is satisfied. However, we point out 
that the bandwidth of the defect level in the A2 and B2 
cases is very sensitive to the interlayer distance, since the 
H atom lies between two adjacent sheets. To check this 
effect, we performed test calculations with variations of 
±5% in this distance and we found that a spin-split can 
also be observed in these cases, in a similar fashion to 
the ABC trilayer, with fractional magnetizations of up 
to O.SjiiB [36]. Since the PBE functional we used in this 
work does not include Van Der Waals interactions, the 
relaxed value for the interlayer distance may be different 
from the one reported here and this effect may be im¬ 
portant. On the other hand, test calculations in other 
configurations (including the case of vacancies, which we 
will discuss later) indicate that small variations in the 
interlayer distance do not induce any important modifi¬ 
cations in the magnetizations or band structures, so the 
Van Der Waals interactions play no major role in these 
cases. 

The magnetizations observed in Table [T| come from the 
full occupation of the spin up defect level in the Al and 
B1 cases, while no magnetization is observed in the A2 
and B2 cases due to lack of spin polarization (for the 
relaxed interlayer distance). The formation of the C- 
H bond during the adsorption effectively removes one 
electron from the honeycomb lattice, leaving an unpaired 
Pz electron in it. This electron then occupies one of the 
spin-split defect levels and induces the magnetization. 
For this reason, we expect that the corresponding spin 


density should be localized around the H atom. We will 
see in the next subsection that this is indeed the case. 
On the other hand, in the A2 and B2 cases, both spin up 
and down defect levels are equally, but not fully occupied. 
The resulting effect then is a shift of the Fermi energy, 
with a small occupation of conduction band states. 



FIG. 3: Band structures without spin polarization for hydro¬ 
gen adsorption on top of Al and B1 sites in the ABA trilayer. 
The defect levels are much flatter in these cases, thus favoring 
spin polarization, as can be seen in Fig. 

Now we move our discussion to the ABC trilayer. In 
this case, the pure band structure near the Fermi energy 
consists of a pair of bands with cubic dispersion, followed 
by two pairs of higher energy quadratic bands, which are 
degenerate at the F point. The presence of inversion sym¬ 
metry in this trilayer leads to a touching point between 
the valence and conduction bands, which lies outside of 
the F point due to trigonal distortion. The hydrogen ad¬ 
sorption lifts this degeneracy and a small gap (< 0.1 eV) 
is opened between the cubic bands. In all cases, the de¬ 
fect levels exhibit spin polarization and, in the Al and B1 
cases, a small spin polarization is also seen in the cubic 
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bands in the vicinity of the F point. This effect is related 
to the high density of states associated with a cubic dis¬ 
persion, due to the almost flat portion near this point. 
When the Fermi energy lies close to the flat portion, as is 
the case for the A1 and B1 sites, spin polarization occurs. 
Similar to the ABA trilayer, a large spin polarization is 
observed for the defect level in the B1 case, which leads 
to a full magnetization. In the other cases, the spin up 
and down levels almost touch at the F point and a small 
occupation is observed for the spin down level, thus re¬ 
ducing the magnetization. In the A2 and B2 cases, there 
is also a shift in the Fermi level and a small occupation 
of conduction band states, again in similarity with the 
ABA trilayer. 


Finally, to finish this subsection, we now discuss the or¬ 
bital nature of the defect level. To this end, we have per¬ 
formed projected density of states calculations (PDOS) 
in the vicinity of the Fermi energy. In Fig. we show 
the PDOS for one particular configuration: adsorption 
on top of an A1 site in the ABA trilayer. All configura¬ 
tions show similar features, including those of the ABC 
trilayer, so we concentrate our discussion in this case. 
As can be seen on the figure, the main contribution to 
the total DOS (black line) near the Fermi energy comes 
from the Pz orbitals of carbon (blue dash-dotted line), as 
expected. The central peaks are indications of the de¬ 
fect levels, since they are mostly flat and thus strongly 
contribute to the DOS. Moreover, the spin up and down 
peaks are split in energy, indicating spin polarization. 
Therefore, the defect level features observed in Fig. 
are correctly reproduced by the PDOS, which also shows 
that this level has mainly a Pz character, with a small 
contribution from the s orbital of the H atom. 



Energy (eV) 

FIG. 4: Projected Density of States (PDOS) for hydrogen ad¬ 
sorption on top of an A1 site in ABA trilayer graphene. Other 
configurations show similar features. Black, red dashed, green 
dotted and blue dash-dotted lines represent the total DOS, 
contributions from the s orbital of hydrogen, a and pz orbitals 
of carbon, respectively. Positive (negative) values represent 
the spin up (down) PDOS. The Fermi energy is set to zero 
and a gaussian broadening of 0.2 eV was used on the energy 
levels. Contributions from polarized orbitals are not shown. 


C. Spin density and local density of states 


In Fig. [^we show the spin density profiles for selected 
defect configurations, where only the layer containing the 
defect is shown. The A2 and B2 cases in the ABA tri¬ 
layer are not shown in the figure, since they are found to 
be non-magnetic, as already discussed in the last section. 
As we can see, in all cases the spin density has a triangu¬ 
lar pattern with opposite values for each sublattice, and 
it is concentrated near the defect site, in agreement with 
previous calculations on monolayer graphene 0128]. The 
strongest magnitude is observed near the carbon atoms 
which are first neighbors of the defect site, while the mag¬ 
nitude in the defect site itself is small, with the lowest 
value observed in the ABC-A2 case. The ABC-Al case 
shows a particularly concentrated spin density, while in 
the other cases the density spreads further into the su¬ 
percell. 
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when the defect site lies on the middle layer, the LDOS is 
more concentrated near the defect center and the peaks 
do not follow a quasi-localized behavior. On the other 
hand, when it is on an external layer, the LDOS spreads 
more throughout the sheet and, particularly, the ABA- 
A1 case shows three peaks at r < 6 A of roughly the 
same size, which are compatible with a quasi-localized 
behavior. In the ABC-Al case, however, we see from the 
profile that the magnitude of the LDOS at the defect site 
is reduced when compared to the previous case (this can 
also be seen in the spin density profiles of Fig. and 
the graph shows that the first two peaks have slightly 
different magnitudes, so they are not fully compatible 
with a quasi-localized behavior. In all cases, for r > 6 
A (roughly half the supercell dimension), the observed 
peaks could contain contributions from periodic images 
of the defect, so we do not consider them in our dis¬ 
cussion. We also didn’t see any LDOS peaks associated 
with second neighbors, which agrees with the prediction 
that the wavefunction amplitude in sites belonging to the 
same sublattice as the defect site in the monolayer should 
be small (it would be zero if electron-hole symmetry was 
present). 


FIG. 5: Spin density for selected configurations. Red (blue) 
regions indicate spin up (down) excess. The A1 and A2 cases 
share similar features with the B1 and B2 cases in the corre¬ 
sponding trilayers. For clarity, only the layer containing the 
adsorbed H atom is shown in each case. 


Fig. shows the local density of states (LDOS) for the 
A1 and A2 configurations in the ABA (top) and ABC 
(bottom) trilayers. As in the spin density, the results for 
the B1 and B2 sites are very similar to those of the Al 
and A2 sites, so we do not show them here. For each 
case, the LDOS was integrated in energy between the 
minimum and maximum values of the spin up defect lev¬ 
els shown in Fig. so it shows essentially the LDOS 
associated to the localized defect state. In the graph¬ 
ics panels, we plot r^x LDOS versus r, where r is the 
distance to the defect center and the LDOS is averaged 
over the polar coordinate. For a quasi-localized state, 
the associated wavefunction decays as 1/r, as predicted 
by tight-binding calculations and observed in STM 
experiments for vacancies in graphite HO]. Therefore, the 
magnitude of the LDOS peaks in a quasi-localized state 
would decay as 1/r^ and r^x LDOS should be a con¬ 
stant at these peaks. We look for this possibility in our 
calculations. For both stackings, we can see then that. 
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FIG. 6: Localized density of states for the A1 and A2 con¬ 
figurations in the ABA (top) and ABC (bottom) trilayers, 
integrated along the bandwidth of the corresponding spin up 
defect level. For each case, the top panels show the profiles in 
the layer containing the defects and the bottom panels show 
the corresponding plots of r^x LDOS as a function of r, the 
distance from the defect center, where the LDOS is averaged 
over the polar coordinate. 


D. Dependence with the external field for the 
ABC trilayer 

To finish our discussion on hydrogen adsorption, we 
now discuss the effect of the application of an external 
electrical field perpendicular to the layers in the ABC 
trilayer. As we have mentioned before, such field opens 
a tunable gap in the pure trilayer. In the first column of 
Fig.[7l we show the evolution of this gap in the pristine 
supercell for four field values, from top to bottom: E = 0 
(same as in Fig. [^, 0.1, 0.5 and 1.0 V/A. In the other 
columns, we show the corresponding band structures for 
different defect configurations. For small field values such 
as in the second row, we can see that the bands coming 
from the pure trilayer evolve in a way similar as they 
do in the pure supercell, with the opening or increase 
of small gaps. The only exception is a linear-like band 


in the B2 case, which crosses the gap between the cubic 
bands and it is already present in zero field. The defect 
bands remain close to the Fermi energy and they are 
slightly changed near the F point, thus indicating a small 
increase in the hybridization with band levels. These 
changes combined with small shifts in their positions lead 
to small magnetization changes for the A2, B1 and B2 
cases, while the A1 case becomes non-magnetic, as shown 
in Fig. [8| 

For higher field values, the hybridization between the 
defect level and the bands becomes much stronger and 
the connectivity of the bands is modified, also making all 
configurations become non-magnetic. For E = 1.0 V/A, 
we can see remains of the flat dispersion of the defect 
level mixed below the Fermi energy for the A1 and B1 
cases and slightly above it for the A2 case. The change 
in connectivity also leads to the formation of a nearly 
flat midgap state in the B2 case, which is degenerate at 
the F point. This could indicate the existence of a truly 
localized state inside the gap, in the sense that the as¬ 
sociated wavefunction becomes exponentially localized, 
as revealed by previous tight-binding calculations m 
[37]. However, we don’t see any qualitative changes in 
the LDOS, with the pattern remaining equal to the one 
observed in Fig. [^(bottom right). We also point out that 
the A1 and A2 cases also show midgap-like states, but 
they have large bandwidths and are very different from 
the original defect levels, again due to rehybridization. 
Finally, we remark again that the difference between the 
A2 and B2 cases, sites that are originally equivalent by 
symmetry in the pure trilayer, lies in our choice of place¬ 
ment of the hydrogen atom. In the former, it lies be¬ 
tween two carbon atoms in adjacent sheets (the B1 and 
A2 sites), while in the latter it lies between the B2 site 
and a hexagon center in the top sheet, thus leading to 
different calculated band structures. 



E (V/A) 


FIG. 8: Evolution of the supercell magnetization with elec¬ 
trical field. Black, red, green and blue lines represent the Al, 
A2, B1 and B2 defect configurations, respectively. 

















FIG. 7: Evolution of the band structure with electrical field for different H adsorption configurations on the ABC trilayer. 
From top to bottom, the rows correspond to field values E = 0, 0.1, 0.5 and 1.0 V/A. The columns correspond to the defect 
configurations indicated in the first row. Black (red dashed) lines represent spin up (down) bands and the Fermi energy is set 
to zero in all cases. 


IV. ISOLATED VACANCY 

A. Structural properties, formation energies and 
magnetic moments 

Now we move our discussion to the case of an isolated 
vacancy. Since many of the results we show here are sim¬ 
ilar to the case of hydrogen adsorption, we will focus our 
discussion on the differences between the two cases. We 
begin with the structural properties, which already show 
important differences. In Fig. we show the structure 
of a layer containing the vacancy after the relaxation. 
The results are very similar for all defect configurations 
and both stackings. The removal of a carbon atom from 
the lattice leaves three unpaired a electrons, one in each 
neighboring atom (labeled 1, 2 and 3 in the figure), along 
with an unpaired Pz electron already seen in H adsorp¬ 
tion. This leads to a larger structural deformation, where 
the atoms 2 and 3 move closer to each other and form 
a weak bond, thus saturating their unpaired electrons. 


The atom 1, however, remains with an unpaired a elec¬ 
tron which, together with the Pz electron, contributes to 
the magnetization. The local environment of the defect 
becomes an isosceles triangle, thus breaking the triangu¬ 
lar symmetry of the lattice. Such a distortion is known 
as a Jahn-Teller distortion, and it has already been iden¬ 
tified in single layer graphene. The bond lengths for spin- 
polarized calculations are reported in Table [TT| and are in 
good agreement with previous calculations on the mono- 
layer, as well as the other properties reported BMM- 
l32] . Here we have not considered the B2 case in the 
ABC trilayer, since the A2 and B2 sites are equivalent 
by symmetry in the pristine trilayer and as such we do 
not expect any noticeable differences between them. For 
unpolarized calculations only, we have also observed an 
out-of-plane distortion in atom 1, which moves about 0.5 
A in that direction, with the exception of the A2 and B2 
cases in the ABA trilayer, where the distortion remains 
purely planar. Such a difference between spin polarized 
and unpolarized calculations has also been reported pre- 
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viously j9l[29]. 

Of course, there are three equivalent distorted con¬ 
figurations, since the atoms 1 and 2, or 1 and 3 could 
also form the bond, leaving the remaining atom with 
the unpaired electron. This could lead to a dynamical 
Jahn-Teller effect, where the system could be constantly 
changing between these configurations, if the energy bar¬ 
rier between them is low enough compared with the tem¬ 
perature [33]. 



FIG. 9: Layer containing the vacancy after structural relax¬ 
ation. The atoms 2 and 3 move closer to each other and form 
a weak a bond, while the atom 1 remains with an unpaired a 
electron, which also contributes to the magnetization. 


TABLE 11: Results for spin-polarized calculations of an iso¬ 
lated vacancy on trilayer graphene: bond lengths of the first 
neighbors of the defect site (atoms 1, 2 and 3 in Fig. |^, forma¬ 
tion energies {Ef) and magnetic moments per unit supercell 
(M) for each configuration studied. For all cases, di 2 = dis. 


Defect Site 

^23 (A) 

di2 (A) 

Ef (eV) 

M (Ats) 

ABA trilayer 

A1 

1.95 

2.56 

7.80 

1.46 

A2 

1.91 

2.54 

7.81 

1.34 

B1 

1.93 

2.55 

7.83 

1.43 

B2 

1.94 

2.55 

7.77 

1.43 

ABC trilayer 

A1 

1.98 

2.55 

7.81 

1.56 

A2 

1.90 

2.54 

7.80 

1.41 

B1 

1.92 

2.54 

7.84 

1.42 


In Table |TT| we have also included the formation ener¬ 
gies and supercell magnetizations for each vacancy con¬ 
figuration. In this case, the formation energies are de¬ 
fined as 

N -I 

^trilayer-\-vac ^trilayer-) (2) 

where Etriiayer+vac Is the total energy of the supercell 
containing the defect, Etriiayer is the total energy of the 
pristine supercell and N is the number of atoms in that 
supercell (in our case. A/" = 72 in a 6 x 6 supercell). 
Comparing Tables |l| and [ll] we see that the vacancies 


have a much smaller sensitivity to the site where they are 
placed, with values that differ by a few tenths of meV. 
This could be related to the large structural distortion 
they induce, which modifies the local environment in the 
layer and is also insensitive to the defect site. However, 
for unpolarized calculations, the formation energies are 
0.25 — 0.5 eV higher and show a larger layer sensitivity, 
which is related to the presence (or lack) of out-of-plane 
distortion in atom 1, as discussed above. 

On the other hand, the magnetizations show a large 
sensitivity to the defect site. Sites that have one C atom 
right on top of them (B1 in both trilayers and also the A2 
site in the ABC trilayer) have roughly the same magne¬ 
tization, while the A1 sites, which have hexagon centers 
on top of them have larger magnetizations. The B2 site 
in the ABA trilayer, which has one hexagon center above 
and one below it has a magnetization similar to those of 
the first group, while the A2 site of the same trilayer, 
which has one site above and one below it, has the small¬ 
est magnetization. This indicates that the mechanism 
responsible for the observed trends could be a charge 
transfer between the layers, which changes the occupa¬ 
tion of the defect levels and reduces the magnetization. 
The interlayer coupling is stronger in sites that are on 
top of each other, so the charge transfer could be more 
intense in these cases (and specially intense in the ABA- 
A2 case). We now discuss these results in terms of the 
band structure. 


B. Band structure 

In Fig. we show the band structures for each va¬ 
cancy configuration. As in the hydrogen adsorption case, 
we can see the presence of defect levels near the Fermi 
energy, which are mainly of Pz character, as we shall see 
below. We can also see the presence of a flat spin up 
level below the Fermi energy, between —0.8 and —0.7 
eV, which was absent in the adsorption case. This level 
has (7 character, and comes from the unpaired a electron. 
It suffers a large exchange split, with the corresponding 
spin down level lying more than 1.0 eV above the Fermi 
level. Therefore, this spin up a level contributes a full 
IpB to the magnetization. The other two unpaired a 
electrons also yield defect levels but, since they form a 
weak bond, these levels he far apart from the Fermi level, 
the exchange splittings are small and hence they do not 
contribute to the magnetization. The remaining contri¬ 
bution to the magnetization comes from the Pz levels, 
which have larger bandwidths than those observed in the 
adsorption case. This indicates a larger interaction with 
band states, leading to charge transfer between band and 
defect levels, partial occupation of both spin up and down 
levels and a reduced contribution to the magnetization 
(< I.Opb) in all cases. In particular, the different oc¬ 
cupations of the spin down pz level in each case are the 
main responsible for the trends observed in Table [HI This 
level is occupied in expense of both spin up Pz level and 
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valence band states, leading to a Fermi level shift (with 
respect to the value in the pure trilayers), which in turn 
could indicate a charge transfer between the layers, as we 
have mentioned before. 

We now discuss the behavior of the Pz levels in each 
case in more detail. In the ABA trilayer, we can see 
that when the vacancy lies on the middle layer (A2 or 
B2 sites), the linear bands from the pristine trilayer are 
preserved, while a sizable gap of 0.1 — 0.2 eV is observed 
for the quadratic bands. However, when the vacancy is 
on a external layer, the inverse situation happens, with 
the quadratic bands preserved and a large gap in the 
linear bands. This could be related to the mirror sym¬ 
metry of the C Pz orbitals in each case. Tight-binding 
calculations in the pristine trilayer show that A2 and B2 
orbitals only contribute to the quadratic bands, while A1 
and B1 orbitals contribute both to linear and quadratic 
bands, by combining with A3 and B3 orbitals with odd 
or even combinations with respect to mirror symmetry, 
respectively m- Therefore, a vacancy in the A2 or B2 
sites should not have an important effect on the linear 
bands, while a vacancy in the A1 or B1 sites may not 
affect the shape of the quadratic bands near the Fermi 
level. This result is in contrast with we have seen on the 
H adsorption. There, a defect on the middle layer also 
leads to a widening of the gap between the linear bands 
([^ , which we believe is related to the small sp^ — sp^ re¬ 
hybridization induced by the formation of the C-H bond, 
thus breaking the mirror symmetry in the vicinity of the 
defect. 

In the ABC trilayer, a similar situation happens. The 
presence of a vacancy in the middle layer does not open a 
large gap between the cubic bands, while a 0.1 eV gap can 
be seen in the A1 and B1 cases. Again, this is related to 
the symmetry of this trilayer. In this case, tight-binding 
calculations show that the cubic bands are formed by 
combinations of Pz orbitals coming from the A1 and B3 
sites, where the interlayer coupling is weaker and, as such, 
they are called low-energy sites m- Therefore, the pres¬ 
ence of a defect in the middle layer should not affect the 
low-energy bands. 

For unpolarized calculations, the results are quite dif¬ 
ferent. In Fig. El we show an example for the Al 
site in the ABC trilayer. As we can see, in this case 
the a defect level coming from the unpaired electron in 
atom 1 lies close to the Fermi level. Due to the energy 
proximity to the Pz level, these levels rehybridize and 
form the observed structure, with larger bandwidths than 
those observed in spin-polarized calculations. This anal¬ 
ysis is confirmed by the PDOS calculation shown in Fig. 

where we compare polarized (left) and unpolarized 
^ght) calculations for the same case. In the polarized 
case, the levels he far apart and thus retain their pure a 
and Pz character, while in the unpolarized case the two 
defect levels near the Fermi energy show contributions 
from both orbitals, as can be seen by comparing a (red 
dashed) curves in each case. As we have mentioned be¬ 
fore, in most cases of the unpolarized calculations the 


atom 1 suffers an out-of-plane displacement which, as in 
the H adsorption case, promotes a local sp‘^ — sp^ rehy¬ 
bridization, thus allowing the mixing between a and Pz 
levels. When there is no such displacement, as in the 
A2 case of ABA trilayer (Fig. 13), these levels are de¬ 
coupled and have smaller bandwidths, compatible with 
those observed in the polarized calculation. 



FIG. 11: Unpolarized band structure for a vacancy on an 
Al site in the ABC trilayer. The out-of-plane displacement of 
atom 1 in this case leads to rehybridization between the a and 
Pz levels, resulting in defect levels with increased bandwidths. 
Other cases show similar features, with the exception of A2 
and B2 cases in the ABA trilayer, shown in Fig. 



E (eV) E (eV) 


FIG. 12: PDOS for a vacancy on an Al site in ABC tri¬ 
layer graphene, comparing spin-polarized (left) and unpolar¬ 
ized (right) calculations. Black, red dashed and blue dotted 
lines represent the total DOS and contributions from a and 
Pz orbitals of carbon, respectively. Positive (negative) values 
represent the spin up (down) PDOS. The Fermi energy is set 
to zero. The peak and shoulder structure on the a curve on 
the right indicates mixing between the a and pz defect levels 
in unpolarized calculations. 
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ABA 



FIG. 10: Band structure results for an isolated vacancy on ABA (top) and ABC (bottom) trilayer graphene. The vacancy is 
on the site indicated in each frame. Spin up (down) levels are represented by black (red dashed) lines and the Fermi level is 
set to zero in all cases. Energies are in eV. For reference, the band structures of the pristine trilayers are shown in Fig. |^(hrst 
column). 



FIG. 13: Unpolarized band structure for a vacancy on an A2 
site in the ABA trilayer. The a and pz levels are decoupled in 
this case due to the lack of out-of-plane displacement of atom 
1, and the corresponding bandwidths are smaller. A second 
cr level, coming from the weak bond between atoms 2 and 3, 
can also be seen within the energy window in this case. 


C. Spin density and local density of states 

We now turn our attention to the induced spin-density. 
Our calculations show that the patterns are insensitive to 
the defect site, and they are all equivalent to the pattern 
shown in Fig. for the A1 case in the ABA trilayer. 
This pattern is very different from the triangular pattern 


observed in the H adsorption case. In this case, the spin- 
density is very concentrated near atom 1, which contains 
an unpaired a electron, as we have mentioned before. 
This is an indication that the contribution coming from 
the (7 defect level to the magnetization is stronger than 
the one from the Pz levels, which we would expect to in¬ 
duce a triangular pattern similar to the one seen in H 
adsorption. As we saw in Fig. the a level is prac¬ 
tically insensitive to the defect site, which explains the 
observed insensitivity in the spin-density. This pattern 
is similar to the one observed for vacancies in single layer 
graphene [iEiiso]. 



FIG. 14: Spin density for the A1 case in the ABA trilayer. 
All other conhgurations show the same pattern. Red (blue) 
regions indicate spin up (down) excess. For clarity, only the 
layer containing the vacancy is shown. 
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On the other hand, the local density of states for the Pz 
levels in spin-polarized calculations shows a strong layer 
sensitivity and is roughly stacking independent. In Fig. 


15, we compare the results for the A1 and A2 cases in the 


ABA trilayer, which are similar to the B1 and B2 cases, 
respectively. The LDOS is integrated along the band¬ 
width of the Pz-^P defect level, which also contains the 
down level within the same energy range (note that the 
inverse situation occurs in H adsorption, as can be seen 
in Fig. so only the spin-up level was included in that 
case). When the defect is on an external layer, the LDOS 
shows a triangular pattern, which is more spread out in 
the supercell than the patterns observed in H adsorp¬ 
tion (see Fig. [^. The corresponding plot of r^x LDOS 
also shows more structure and it is roughly constant for 
1.5 A< T’ 4 A, which could indicate a Quasi-localized 
behavior in this region. For intermediary distances this 
behavior is not well satisfied and, for larger distances, 
larger peaks are observed and they may contain contri¬ 
butions from the periodic images, as was also pointed out 
in the H adsorption case. 





FIG. 15: LDOS for the A1 and A2 configurations in the 
ABA trilayer, integrated along the bandwidth of the pz defect 
levels (both spin up and down). For each case, the top panels 
show the profiles in the layer containing the vacancy and the 
bottom panels show plots of r^x LDOS as a function of r, in 
the same fashion as Fig. 


When the defect is in the middle layer, the situation is 
very different. The LDOS show s a p attern that resembles 
the spin-density shown in Fig. M The r^x LDOS plot 
shows a two-peaked structure around atom 1, and a small 
peak around the bond between atoms 2 and 3 can also 
be seen in the profile map. Therefore, the defect level 
wavefunction is more localized near the vacancy center 
and does not follow the 1/r^ power-law. We have seen a 
similar situation in the hydrogen adsorption case, where 
the LDOS is also more localized at the defect center when 
it is on the middle layer, although the observed patterns 
are very different. Note also that the LDOS amplitudes 
of the spin up and down levels (black and red lines in 
the line profiles) are very different near the peak at atom 


1, but they are very similar for larger distances. This 
could indicate that this long range contribution, which 
decays weaker than 1/r^ (since r^x LDOS increases with 
r) and has a very weak magnitude when compared to the 
previous plots, could come from band levels within the 
integration window, which would add to possible contri¬ 
butions from periodic images. In the other plots, this 
contribution may be present, but its magnitude is much 
weaker than the others we have discussed previously. 


Finally, we mention that, for unpolarized calculations, 
we do not see a layer sensitivity and all patterns are sim¬ 
ilar to the localized pattern we have just discussed. This 
could be related to the presence of a cr level near the 
Fermi energy in this case which, as we have seen, rehy¬ 
bridizes with the Pz level when the defect is on an exter¬ 
nal layer. The corresponding wavefunctions then have a 
contribution from the a orbital coming from the unpaired 
electron, leading to a stronger localization. 


D. Dependence with the external field for the 
ABC trilayer 


We now discuss the effect of the application of an ex¬ 
ternal electrical field in the ABC trilayer with a vacancy. 
In Fig. we show the band structures for zero field 
and three field values for each vacancy configuration, as 
we did in Fig. The trends observed here are very 
similar to those of H adsorption. For a small field value 
(second row), the band levels remain similar to those of 
the pure trilayer, but the gap between the cubic bands 
is affected by the vacancy. For the A1 and B1 configu¬ 
rations, the gap is reduced and almost closes, indicating 
that the charge transfer between the layers induced by 
the field is compensated by an opposite charge transfer 
due to presence of the vacancy (the field is applied in the 
direction of layer 1 to layer 3). A similar behavior is seen 
on bilayer graphene, when potassium is adsorbed on top 
of a layer when the system has residual doping [34]. For 
the A2 case, however, the gap is very small at zero field 
and it evolves naturally as the field increases. The small 
gap at zero field in this case indicates that the external 
layers contribute almost equally to the charge transfer to 
the middle layer in this case, so they have roughly the 
same electrical potential. When they have the same po¬ 
tential, no gap is induced, in agreement with our earlier 
discussion in section IV B in terms of the symmetry of 
the ABC trilayer and the Pz orbitals. 







13 



FIG. 16: Evolution of the band structure with electrical held for different vacancy conhgurations on the ABC trilayer. From top 
to bottom, the rows correspond to field values E = 0, 0.1, 0.5 and 1.0 V/A. The columns correspond to the defect configurations 
indicated in the first row. Black (red dashed) lines represent spin up (down) bands and the Fermi energy is set to zero in all 
cases. Band structures for the pure trilayer are shown in Fig. 



For larger field values (third and fourth rows), the Pz 
defect levels strongly mix with the band levels, changing 
the band connectivity and reducing the magnetization 
(Fig. 17). For F; = 0.5 V/A, we can see a ’’midgap” 
state with a large bandwidth in the Al and B1 cases 
and a small one in the A2 case, but none of them lie 
completely inside the gap between the adjacent bands. 
By increasing even more the field, the spin polarization 
eventually disappears and the sigma level is pulled back 
to the Fermi level, as we have discussed in section EXl 
and can be seen in the last row. The Pz defect level 
(now rehybridized with band levels) remains visible only 
in the A2 case, and still has a small bandwidth. The 
corresponding LDOS for this level has the same localized 


FIG. 17: Evolution of the supercell magnetization with elec¬ 
trical field. Black, red and blue fines represent the Al, A2, 
B1 defect configurations, respectively. 
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pattern shown in Fig. (right), so this case could be a 
candidate for a truly localized midgap state, as predicted 
by a previous tight-binding calculation HZ]. We remark 
that a similar situation happened in H adsorption for the 
B2 case, as we have discussed in section alD[ We expect 
this case to actually resemble more the result of a typical 
tight-binding calculation since, as we have discussed, H 
adsorption induces small structural changes and remove 
only a Pz electron from the lattice (tight-binding calcu¬ 
lations in graphene and related systems usually neglect 
structural distortions and include only Pz electrons). 

V. CONCLUSIONS 

We have studied the electronic and structural proper¬ 
ties of hydrogen adsorption and single vacancies on tri¬ 
layer graphene. We have verified that most properties, 
such as formation energies, magnetization and local den¬ 
sity of states, show a sensitivity to the layer in which the 
defect is placed, followed by smaller sensitivities to sub¬ 
lattice and stacking type. Both types of defects remove a 
Pz electron from the lattice, thus inducing a defect level 
in the band structure near the Fermi energy. This level 
shows a bandwidth, as a result of the interaction with the 
band levels near the F point of the supercell BZ (both K 
and K' points from the primitive cell BZ are folded into 
this point), which is larger for the vacancy case. For spin 
polarized calculations, this level also suffers an exchange 
splitting due to the zero energy instability (large density 
of states at the Fermi energy due to the flat portion of the 
level), and the different populations of the spin up and 
down levels lead to the different magnetizations observed 
in each defect configuration. For vacancies, the removal 
of a carbon atom also leaves three unpaired a electrons. 


two of which form a weak bond, while the third remains 
unpaired, resulting in a Jahn-Teller structural deforma¬ 
tion. These electrons also induce defect levels, one of 
which lies close to the Fermi level (the one coming from 
the unbound electron) and suffers a large exchange split¬ 
ting, thus also contributing to the magnetization. Un¬ 
polarized calculations may show a rehybridization be¬ 
tween this level and the Pz level, depending on whether 
an out-of-plane displacement of the atom containing the 
unpaired electron occurs. 

Spin density calculations exhibit a triangular pattern 
in H adsorption and a localized pattern in vacancies, 
where the localization occurs near the atom containing 
the unpaired a electron, whose contribution is stronger 
than that of the Pz defect level. Local density of states 
show patterns that are more localized near the defect 
center when the defect is placed on the middle layer. Fi¬ 
nally, our calculations with an applied external electrical 
field in the ABC trilayer indicate that the Pz defect level 
strongly mixes with band levels for high field values, for 
both type of defects, resulting in reduction and eventu¬ 
ally loss of magnetization. There is also a possibility of a 
midgap-like state with a small bandwidth when the defect 
is placed on the middle layer, as predicted by previous 
tight-binding calculations. Such a state can be very im¬ 
portant for the transport properties of the ABC trilayer 
when used in electronic devices. 
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